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Abstract 

Binding energies and density profiles of two-dimensional systems of liquid 4 He with different 
geometries are studied by means of a zero-range density functional adjusted to reproduce the line 
tension obtained in a previous diffusion Monte Carlo calculation (Xdmc = 0.121 K/A). It is shown 
that this density functional provides accurate results for the binding energy of large clusters with 
a reasonable computational effort. 

PACS numbers: 31.15.Ew 36.40.-c 61.46. +w 67.70.+n 68.03.Cd 
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I. INTRODUCTION 



Quantum liquids in restricted geometries have attracted a lot of attention in recent 
years. One interesting feature of these systems is that their internal structure is more 
accesible than in bulk liquids due to the restricted motion of the atoms in the confining 
potential. Among these systems the study of two-dimensional systems has received partic- 
ular attention. An example of such a system is liquid helium adsorbed to a more-or-less 
attractive flat surface. This kind of system was observed for the first time by M. Brctz et 
al. jsj in 1973, when they reported the observation of adsorbed 4 He onto the basal plane of 
graphite. In the last few years, adsorption properties of helium on many different substrates 
(carbon, alkali and alkaline-earth flat surfaces, carbon nanotubes, aerogels) have become a 
fruitful topic of research. 

Theoretical microscopic studies with realistic atom-atom interactions, such as that of 

n 

Clements et al. [3] , have shown that films with low surface coverages, where all atoms cover 
the surface with a thickness corresponding to a single atom, can be reasonably approximated 
by a two-dimensional model. In connection with these systems, an interesting question 
naturally arises as how physics depends on the dimensionality of the space. 

The homogeneous 2D liquid has been studied using different theoretical methods, such as 
molecular dynamics j^J and quantum Monte Carlo simulations either Green's Function 
or diffusion p techniques. 

Recently, two-dimensional clusters of liquid 4 He have been studied using using a shadow 
variational wave function [ij, and also by the diffusion Monte Carlo (DMC) method 
In these two references, the binding energies of the 2D-clusters were fitted by means of a 
mass formula and a line tension of A = 0.121 was reported in Ref. 8j]. However, due to 
computational limitations, the number of atoms in the clusters was limited to N ~ 100. 
Besides, the density profiles of the clusters, specially for coordinates close to the origin, 
are usually obtained in the Monte Carlo method with poor statistics. Therefore it seems 
appropiate to build a density functional suitable to be used in 2D, using the same procedure 
which has already been succesfully used to study 3D 4 He clusters joL . 

Density functionals are based on the well-known Hohenberg-Kohn theorem [ljj that as- 
serts that the ground-state energy per particle of a many-body system can be written as 
a functional of the density. Once the functional is available, its minimization brings to an 
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Euler-Lagrange equation for the density profile p(r), which allows to calculate the properties 
of clusters. Note that the results obtained with the density functional will be more reliable 
for larger clusters, as it has been constructed to reproduce properties of the homogeneous 
and semi-infinite media. 

The density functional we use is the simplest version of the zero-range functional inten- 
sively used in 3D Its parameters have been adjusted so as to reproduce some properties 
of the ground state of the homogeneous system as obtained in DMC calculations p , as well 
as the line tension extracted from the mass formula of Ref. ^]. This procedure is discussed 
in Sect. HH together with the results for the slabs. Sect. IIHI is devoted to the study of 
finite droplets, with special emphasis for those with a large number of atoms. Finally, in 
Section HVI the main conclusions are summarized. 



II. SEMI-INFINITE SYSTEM AND SLABS 

Density functionals to investigate surface properties of superfluid 4 He were developed 
during the 1970's [14]. At zero temperature and in the absence of currents, the order 
parameter of a bosonic system is nothing but the square root of the one-body density 
p(r), so it is natural to think of the energy of the system as a functional of the helium 
density. Stringari ^| proposed a zero-range density functional for non- homogeneous three- 
dimensional 4 He systems of the form: 

E[p] = J dr + hp- + cf* + rf|V P | 2 } , (1) 

which was used to study 4 He surface properties [3] and clusters [^J. 

For the two-dimensional systems, we will use the same functional form. However, the 
parameters b, c, 7 characterizing the functional for homogeneous systems have to be re- 
calculated by requiring that the functional reproduces some known properties of the two- 
dimensional 4 He system, such as the ground-state energy per particle (eo = —0.89706 K), 
saturation density (po = 0.04344 A -2 ), and the compressibility (or the velocity of sound, 
s = 92.8 m/s), which were derived in the framework of the Difusion Monte Carlo method jf|. 
In addition, the parameter d is obtained by demanding that the line tension of the semiinfi- 
nite system equals that obtained from DMC calculations for 2D clusters, A = 0.121 K/A 
The energy per particle for the homogeneous system provided by this functional reproduces 
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very well the equation of state obtained from DMC calcualtions in a wide range of densi- 
ties H. Therefore, as a first step we need to study the semi- infinite system. In principle, 
one should solve an Euler-Lagrange equation for the density profile p(r) which results from 
minimizing the energy functional of Eq. ([!} after introducing the Lagrange multiplier (p), 
which is identified with the chemical potential, 
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One should also impose the boundary condition, p(x — > — oo) = po- However, for the 
particular case of a zero range functional, the line tension of the semi-infinite system can be 
evaluated in a closed form, without solving for the density profile: 

po r/*2 \ ii/2 
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+ pdj (bp + cp 1 ^ - p) 

Then, imposing A = 0.121 K/A H one gets an implicit equation A = \(d) which can be 
solved numerically, thus fixing the last parameter of the functional. The parameters for the 
two-dimensional zero range functional are listed in Table |U 

Once the functional is defined, one can study any two-dimensional 4 He system. The first 
systems we have considered are two-dimensional slabs with varying central density p c = p(0). 
A systematic study of three-dimensional 4 He slabs with different density functionals has been 
presented in Ref. 16[ The Euler equation for the slabs is the same as for the semi-infinite 
system (Eq. (J2J)). The changes in the solution of the equation originate from the different 
geometry and boundary conditions which define the slab. Translational symmetry implies 
that the density depends only on the coordinate perpendicular to the slab surface, which we 
call x. Then, Vp = p'(x) and V 2 p = p"(x), where the prime denotes derivative with respect 
to x. In this way, the Euler equation (Eq. (J2J) can be expressed in a more convenient form, 
in which now p depends only on x, 
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In the next step, one eliminates the second derivative of Eq. (jlj) by multiplying both sides 
of the equation by p' and integrating with respect to x, from the origin to a given value of 

x, 
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Imposing p(oo) = p'(oo) = 0, and considering that the slab is symmetric respect to x = 0, 
and therefore p'(0) = 0, one obtains the chemical potential as a function of the central 
density of the slab, 

p = bp c + cp< +1 . (6) 

We remark that this chemical potential is constant along the profile. 
Going back to Eq. (0), and using the fact that p'(0) = 0, one obtains 



f ?(bp + cprr+ 1 - f i) 
Y h 2 /8m + dp ' 1 ' 

This expression for p 1 is then used to calculate the number of atoms per unit length along 

the y axis (parallel to the surface), i.e. the coverage, in a closed expression in terms of p, 
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This integral has a singularity when p — > p c , which can be avoided performing and integration 
by parts. The final expression, free of numerical problems reads 
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Also useful is the energy per unit length, 
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Using the Euler equation (Eq. EJ) and the previous definition of the coverage, one can finally 
express the energy per particle e = E/N in terms of the inverse of the coverage, x = L/N, 



r (0) u n 2 \ 

e(x) = p{x) +4x J dp J f — + dp] {bp + cpi+T - p(x)) (10) 

Therefore, given p(0) G [0,p ], and using the previous expressions, one can calculate the 
chemical potential, the coverage and the energy per particle. 

The energy per particle and the chemical potential as a function of the inverse of the 
coverage x, are reported in Fig. (0). In the limit x — > one recovers the binding energy at po 
of the uniform system, which in turn coincides with the chemical potential. The energy per 



particle has a very clean linear behavior at the origin as illustrated by the solid straight line 
which provides a very good description of the energy per particle up to x ~ 1.5 A. The slope 
of this line can be analytically derived and turns out to be twice the linear tension. As a 
consequence of the linear tension, the binding energy per particle of the slab decreases with 
the inverse of the coverage. Actually, the derivation of the linear behavior of the binding 
energy per particle can be obtained starting from Eq. (|1U|). which defines the energy per 
particle as a function of x, by performing an expansion around 5 = 0, 

e(x) = Hoc + 2ax + ... (11) 

For values larger than x ~ 1.5 A the binding energy per particle starts to bend towards 
the x axis and becomes a convex function, which will slowly approach zero. The chemical 
potential is very flat at the origin being determined by the central density of the slab. 
The ratio of the central density to the equilibrium density as a function of the inverse of 
the coverage is displayed in the lower panel of Fig. (0). In agreement with the chemical 
potential, the central density is very flat at small values of x. The central density of a slab 
can never go above the equilibrium density and it is always a decreasing function of x. The 
flatness of the central density and the chemical potential for small values of x, indicates that 
the slab approaches very slowly the limit of the infinite system. 

The density profiles of the slabs can be obtained from the following relation 
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valid for x > 0. Notice also the presence of a divergency when p — > p c which can be 
again avoided performing an integration by parts. The profiles calculated for various central 
densities are plotted in Fig. (J2J). The size of the slab increases with the central density. A 
measure of this size is given by the radius R, defined as the distance from the origin to the 
point where the density has fallen to half its central value. The radius of the slabs as a 
function of x is shown in the top panel of Fig. @. As expected, the radius diverges in the 
limit p c — > po and is a decreasing function of x. It also presents a very shallow minimum 
around x ~ 3.3 A. The profiles are also characterized by the thickness t, defined as the 
distance between the points where the density has decreased from 90% to 10% of its central 
value. The thickness as a function of x is reported in the lower part of Fig. (jHJ). Up to 
x ~ 1, the thickness is a flat function, indicating that the surface of the slab is very much 



the same, actually is the radius of the slab that grows very fast and diverges when x — > 0. 
The thickness is very large for the smaller slabs: as each particle interacts with very few 
others, the system extends to very large distances. When the central density is increased (x 
decreases), the thickness decreases until it has a minimum at x « 2.75 A, with t » 8.6 A. 
Then, it increases again approaching a finite value corresponding to the semiinfinite medium, 
t « 11.03 A. 



III. DROPS AND LINE TENSION 



As a next step, we consider finite systems and, in particular, drops of fixed number of 



atoms N. These where already studied by DMC techniques in Ref. |8(. However, compu- 
tational limitations allowed only to study small values of N. Here, we will take advantage 
of the computational feasibility of the Density Functional calculations and will extend the 
analysis to much greater N. In this way it is possible to study the asymptotic behavior of 
several quantities which characterize the drops. 

In this case, the equation (J2J) for the profile can be rewriten in the form of a Schrodinger- 
like equation for p: 
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As the number of atoms is a well defined N, one would need to be very careful in 
determining the central density p(0) so that the chemical potential adjusted exactly to N. 
However, it turns out that this equation is more efficiently solved by means of the steepest 
descent method An initial trial p(r) is projected onto the mininum of the functional by 
propagating it in imaginary time. In practice, one chooses a small time step At and iterates 
the equation 

p(r, t) « p(r, t) - AtHp(r, t) (14) 

by normalizing p to the total number of atoms at each iteration. The time step that governs 
the rate of convergence should be taken appropriately small in such a way that Eq. (|14j) is 
a valid approximation. Convergence is reached when the chemical potential has a constant 
value independent of position. 

The energy per particle (empty circles) and the chemical potential (full circles) of each 
drop are reported in Fig. (jlj) as a function of N~ l l 2 . Also shown are the DMC results 



(empty squares) and their quadratic fit reported in Ref. |8j. The calculated energies of the 
droplets can be represented very accurately with a mass formula of the type 

e(N) = e b + e lZ + e c z 2 + ..., (15) 

with z = N~ l l 2 . The two first coefficients of this expansion are the bulk energy e& and 
the line energy ej, out of which the line tension A is defined by 27rr A = q. Here r is 
the unit radius, defined as the radius of a disk whose surface is equal to the inverse of the 
equilibrium density of the 2D bulk liquid, i.e, poiir 2 , = 1. Finally, e c can be related to the 
so-called curvature energy. Contrary to the DMC calculations where the largest droplet that 
we studied had 121 atoms, here we have considered droplets with up to 10000 atoms. In 
this way we can accurately study the behavior of the energy per particle for small values of 
z. Doing a quadratic fit to the calculated energies per particle for N > 512 and including 
also the bulk binding energy for z = one gets 

y = -0.897 + 2.0587,2 + 0. 66466z 2 , (16) 

which is plotted by a solid line in Fig. (jlj). One sees that accurately reproduces the bulk 
energy per particle, which was used to fix the parameters of the functional. The value of 
ei = 2.0587 corresponds to a line tension A = 0.1217^/A, which is the same as the value of 
the line tension of the semiinfinite system, used to build the density functional. Note that 
this fit, even if it has been calculated for N > 512, is rather accurate down to N = 36. 
Obviously, one can not expect a good agreement for N = 16. The linear behavior of the 
chemical potential as a function of z is easy to understand using the mass formula (Eq. (|15jl ) 
and the thermodynamic definition of the chemical potential fi = dE/dN, where E is the 
total energy of the droplet. Using this prescription, the slope of the chemical potential as 
a function of z results to be ej/2. A similar plot for the energy per particle in the three 
dimensional case in terms of iV -1 / 3 , would provide a behavior of the chemical potential for 
large N, dominated also by a linear component with a slope at the origin given by 2e s /3, 
where e s is the surface energy associated to three dimensional clusters 9|. 

Also interesting is the fact that the coefficient of z 2 is positive. This signe corresponds 
to the expected loose of binding energy associated to the curvature of the contour of the 
cluster. This is in contrast with the value of e c obtained by fiting DMC results as it was 
done in Ref. [gj . However, in that case the number of particles in the clusters used to build 
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the fit was much smaller, being N = 121 the largest number of particles and going down to 
N=8 for the smallest one. In the present fit we have explicitly avoided the clusters with a 
small number of particles which can easily distort the results, and we have considered only 
the cases with N > 516. 

Next thing to analyze are the density profiles which are reported in Fig. (j5J) for different 
numbers of atoms. Contrary to the slabs, the central density of the droplets can be higher 
than the saturation density, which is indicated in the figure by an horizontal line. The 
profiles are well adjusted by a function of the type: 

p(r) = 7 Pf r ^» ( 17 ) 

1 + e~ 



which has an associated central density p(0) = P//(l + e~ R ^ c ) u . The parameters defining 
the profiles for the different numbers of atoms are provided in Table |HJ together with the 
thickness and root-mean-square radius obtained from these fits. 

The upper part of Fig. © reports the central density of the different droplets as a function 
of z. For large values of N, the central density is larger than the saturation density, i.e. the 
central part of the droplet is more compressed than the bulk system, which is sometimes 
referred to as a leptodermous behaviour |9j. Of course for N — > oo the central density 
tends to the equilibrium density of the homogeneous liquid. First, the central density grows 
almost linearly with z, reaches a maximum, for iV ~ 60, which would correspond to the 
most compressed droplet and then decreases. Finally for N < 25 the central region of the 
droplets becomes less compressed than the bulk system. 

The mean square radius is shown in the lower panel of Fig. (JBJ) as a function of N 1 / 2 . 
The expected linear behavior, associated to a constant average density 

{ rY 2 = -^N^ (18) 

is rather apparent. The fit to the calculated values, from N = 16 to N = 10000, provides 
a p — 0.0434 A~ 2 in very good agreement with the equilibrium density used to define the 
parameters of the density functional. 
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IV. CONCLUSIONS 



We have constructed a density functional suitable to study non-homogeneous two- 
dimensional 4 He systems. First, we have considered two dimensional slabs, to study the 
energy per particle, the central density, radius and density profiles as a function of the cov- 
erage. We have analitically shown that the extracted linear tension from a mass formula 
adapted to this type of geometry, is consistent with the value of the linear tension of the 
semi-infinite system which was used to build the density functional. The thickness and the 
central density of the slabs approach from below the values corresponding to the semi-infinite 
system while the radius diverges. 

We have also studied the energetics and structure of two dimensional clusters. In partic- 
ular, we have considered clusters with a very large number of atoms to study the behavior 
of the mass formula and to establish how the system approaches the bulk limit. The central 
density of the clusters when N — > oo approaches the saturation density from above, and 
therefore the internal regions of the clusters are more compressed than the bulk system while 
the external regions have densities which would correspond to negative pressures or even 
below the spinodal point for a uniform system. 

The profiles of the clusters are very well fitted by a generalized Fermi function. The 
thickness of the cluster slowly approaches the thicknes of the semi-infinite system, as it 
also hapens in the case of the slab geometries, but in this case from above. Finally, we 
have analyzed the linear behavior of the rms radius of the droplets in terms of N 1 ^ 2 , and 
recovered the saturation density from the slope of this fit. The proposed density functional 
can be used with a very small computational effort for large clusters where the value of N 
is prohibitive for a Monte Carlo calculation. 
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TABLE I: Parameters of the two-dimensional zero-range functional. 
b [KA 2 ] c [kA 2 ( 1+ ^] 7 d [KA 4 ] 
-26.35 4.88xl0 5 3.62 359 



TABLE II: Parameters of a generalized Fermi-profile fit (Eq. 117(1 to the density profiles obtained 
with the zero-range density functional. All lengths are in A and pf is in A~ 2 . The parameter v is 
adimensional. 
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pj R c v 
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(r 2 )V2 


16 


0.04321 13.2718 3.22067 2.13417 


11.746 


9.66 


36 


0.04494 19.1852 3.22054 2.37516 


11.548 


12.793 


64 


0.04974 24.6826 3.16845 2.37072 


11.364 


16.28 


121 


0.04441 32.8314 3.12302 2.33636 


11.226 


21.72 


512 


0.04392 64.4245 3.08867 2.32368 


11.112 


43.533 


1024 


0.04378 89.8497 3.08584 2.32997 


11.097 


61.342 


2500 


0.04366 138.628 3.08552 2.33826 


11.090 


95.678 


10000 


0.04355 274.026 3.08653 2.3467 


11.088 


191.276 
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FIG. 1: (top) Energy per particle (dashed line) and chemical potential (dot-dashed line) for 4 He 
slabs as a function of the inverse of the coverage. The solid straight line corresponds to the 
asymptotic behavior of the energy per particle and its slope is determined by the linear tension 
of the semiinfinite system, (bottom) Ratio between the central density and the bulk equilibrium 
density for 4 He slabs as a function of the inverse of the coverage. 
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FIG. 2: Density profiles for 4 He slabs with central densities p c //9 =0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 
0.95, where po is the equilibrium density of the bulk. The horizontal dotted line corresponds to 
the bulk density po- 
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FIG. 4: Energy per particle (empty circles) and chemical potential (full circles) of 4 He droplets as 
a function of N~ l l 2 . Also shown is a quadratic fit (see text) of the results with iV > 516 (solid 
line). The straight short-dashed line is obtained when the mass formula (with terms up to z 2 ) is 
used to calculate the chemical potential. The empty squares are DMC results from Ref. and 
the dotted line corresponds to the quadratic fit to these results reported in the same reference. 
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FIG. 5: Density profiles for 4 He droplets for N=16, 64, 121, 512, 1024, 2500 and 10000 atoms. 
The continuous lines are generalized Fermi profiles (Eq. I17|) fitted to the data (see Table [n]) . The 
dotted horizontal line indicates the equilibrium density pq. 
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FIG. 6: (top) Central density of 4 He droplets as a function of iV -1 / 2 . The empty circles correspond 
to the results obtained with the zero-range density functional, while the line stands for a cubic 
spline fit to these data. The dotted horizontal line indicates the saturation density pq. (bottom) 
Mean square radius of 4 He droplets as a function of N 1 / 2 . The solid line is a linear fit to the data 
without independent term. The empty circles correspond to the results of the zero-range density 
functional. 
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